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ABSTRACT 

We incorporate non-thermal excitation and ionization processes arising from non- 
thermal electrons that result from 7-ray energy deposition, into our radiative trans- 
fer code CMFGEN. The non-thermal electron distribution is obtained by solving the 
Spencer- Fano equation using the procedure of iKozma fe Franssonl (|1992| ). We applied 
the non-thermal ca lculations to the blue sup ergiant explosion model whose early evolu- 
tion was studied in lDessart fe Hillier! (|2010f ). Non-thermal processes generally increase 
excitation and ionization and decrease the temperature of the ejecta. We confirm that 
non-thermal processes are crucial for modeling the nebular spectra. Both optical H 1 
and He 1 lines are significantly strengthened. While optical He 1 lines are not easily dis- 
cerned in observational spectra due to severe blending with other lines, He 1 2.058 /im 
provides an excellent opportunity to infer the influence of non-thermal processes. We 
also discuss the processes controlling the formation of the He 1 lines during the nebular 
epoch. Most lines of other species are only slightly affected. We also show that the in- 
clusion of Fe 1 has substantial line-blanketing effects on the optical spectra. Our model 
spectra and synthetic light curves are compared to the observations of SN 1987A. 
The spectral evolution shows broad agreement with the observations, especially Ha. 
The uncertainties of the non-thermal solver are studied, and are expected to be small. 
With this new addition of non-thermal effects in CMFGEN, we now treat all known 
important processes controlling the radiative transfer of a supernova ejecta, whatever 
the type and the epoch. 

Key words: radiation transfer - radiation mechanisms: non-thermal - stars: atmo- 
spheres - stars: supernovae - stars: evolution 
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1 INTRODUCTION 

Supernovae (SNe) explosions are classified as thermonuclear 
SNe and Core-Collapse SNe (CCSNe), according to their 
explosive mechanisms. Thermonuclear SNe are also called 
Type la SNe, with light curves that can be standardized 
to measure cosmological distances. Based on their spectro- 
scopic signatures, CCSNe can be further categorized into 3 
broad classes, Type lb, Type Ic and Type II. Type II SNe 
exhibit hydrogen lines, while Type I SNe do not. 

SNe explosively produce 56 Ni which decays into 56 Co 
with a half-life of 6.1 days. 56 Co further decays into 56 Fe 
with a half- life of 77 days. Such radioactive decay is the main 
energy source of Type la SNe. While the energy source of 
Type Ib/c SNe resembles that of Type la SNe, Type II SNe 
are first powered by the shock-deposited energy, but at later 
times radioactive decay will be the dominant energy source. 
High-energy 7-rays are produced during the radioactive de- 
cay. The 7-rays undergo Compton scattering with bound 
and free electrons, which results in high-energy electrons. 



These high-energy electrons interact with the medium and 
deposit energy through heating, excitation, and ionization. 

Non-thermal excitation and ionization processes have 
been suggested to explain some observational features of 
various SNe of different types. For example, the Bochum 
event in the Ha profiles of SN 1987A was explained by non- 
thermal excitation due to a n asymmetric ej ection of 56 Ni 
' Phillips fe Heathcotd 1 19891 : IChugail Il99ll ). iLi fe McCra^l 
19951 ) included two-photon processes in models for SN 



1987A, and showed that they could fit the evolution of the 
He 1 1.083 fim and 2.058 /^m lines whose upper levels are pop- 
ulated by non-thermal processes induced by 7-rays. 

Several groups dHarkness et al.| Il987 ; Lucy 

Swartzlll99ll: IKozma fe Franssonl Il992l ; ISwartz et al.l 



1991 



1993; 



Kozma fe Fransson 1998al \h studied Hel lines in Type lb 



SNe a nd suggested the s e line s are due to non-thermal exci- 
tation. lHarkness et al.l (| 1987h carried out LTE calculations 
for Type lb SNe and showed that these did not match ob- 
servations - emission line strengths indicate large departures 
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from LTE for He I. ILucvI (|l99ll ) treated non-thermal excita- 
tion and ionization for He I in their LTE models, and found 
that it could explain the presence and strength of He I lines. 
However, these studies did not treat all atoms in non-LTE 
and/or ignored effec ts due to radiative transfer. Recently, 
IDessart et alJ ()201ll ) carried out a comprehensive study of 
SNe which result from quasi-hydrogen-less low-mass cores 
from 15-25 M main-sequence stars in binary systems using 
CMFGEN. They found that the He I lines could be reproduced 
at early times if the He mass fraction is greater than 50 per 
cent, without invoking non-thermal excitation. 

Line formation in SNe is inherently complex. The fast 
expansion velocities exacerbate the influence of lines, by pro- 
ducing broader lines, increasing line overlap, and enhancing 
line blanketing. Line blanketing strongly influences the op- 
tical and UV ranges of SNe spectra, generally producing a 
flux deficit in those regions. To accurately account for the 
influence of line blanketing requires that an enormous num- 
ber of metal lines be included in the model. Because of non- 
LTE and time dependence, an accurate interpretation of the 
spectral signatures requires quantitative studies. 

For our studies we utilized CMFGEN, which is a radia- 
tive transfer code developed for modeling Wolf-Rayet stars 
(|Hillierl l 19871) that was later e nhanced to allow for line blan- 
keting ([Hillier fc Millerlll998l ). to model al l early- type stars, 
and u ltimately to model all types of SNe (|Hillier fc Dessart] 
I2012D . The code includes the following processes: bound-free; 
bound-bound (including two-photon emission); free- free; di- 
electronic recombination; electron scattering; Rayleigh scat- 
tering by hydrogen; inner-shell ionization by X-rays; col- 
lisional ionization/recombination; collisional excitation/de- 
excitation and charge exchange reactions with H and He. 
The code does full-non-LTE calculations at great computa- 
tional cost. Atomic dat a utilized by CMFGEN are listed by 
IDessart fc Hillierl (|2010l ). 

Time dependent terms in the rate equations are usu- 
ally neglected in SNe spectral modeling at early times, be- 
cause the recombination time is much smaller t han the 
age of the SN. However, lUtrobin fc Chugail (|2005h showed 
that time-dependent ionization of hydrogen is responsi- 
ble for reproducing the strong Ha line s at early times 
in SN 1987A, which was confirmed by IDessart fc Hillierl 
(|2008l ). CMFGEN takes into account the time-dependent 
terms in the statistical and radiative equilibrium equa- 
tions, and in the radiative transfer equations. The code has 
been ext ensively tested and used to s tudy v arious types 
of SNe (IDessart fc Hillierl 120051. [20061. |2008| , |2010| . 1201 ll ; 
IDessart et al.ll2008l ; iHillier fc Dessartll2012l ). 

At early times, non-thermal effects in Type II-P SNe 
are small, since the optical depth is large and the ejecta is 
essentially in LTE in the region where 56 Ni is abundant. 
However, non-thermal effects will eventually take over other 
processe s, especially during the nebular phase. Further, the 
study of Ijerkstrand et ahl ( 201ll ) shows that radiative trans- 
fer effects are still very important at very late times. In 
Type lb SNe, non-thermal processes c a n become important 
as ea rly as 10 days (e.g., ILucvI [l99ll ; IDessart et al.|[201ll . 
120121 ). Thus non-thermal excitation and ionization, as well 
as radiative transfer effects, need to be taken into account in 
order to undertake comprehensive studies of SNe from the 
photospheric through the nebular phase. 

With the implementation of non-thermal processes in 



CMFGEN, we are now able to model from just after the ex- 
plosion, through the photospheric phase, the onset of the 
nebular phase, and the nebular phase using a single code. 
The aim of this paper is to describe the implementation of 
non-thermal excitation and ionization into CMFGEN, and to 
benchmark the non-thermal solver using SN 1987A. The pa- 
per is organized as follows: We present our methodology of 
the non-thermal solver in Section^ and the hydrodynamical 
model, utilized in our studies, is summarized in Section [3] 
Non-thermal effects are discussed in Section[4] while we com- 
pare the properties of the non-thermal and thermal models 
in Section [S] To examine the importance of additional opac- 
ity sources we study the influence of Fel in Section [6] Com- 
parison to the observations is made in Section [7] We explore 
the uncertainties underlying the non-thermal solver in Sec- 
tion^ and discuss the effect of mixing and 7-ray transport in 
Section [5] Finally, we will summarize our results in Section 
1101 A parallel paper will discuss the influence of non -thermal 
processes on Type Ibc spectra |Dessart et al.|[2012l ). 



2 THE METHOD 

Non-thermal excitation and ionization enter the radiative 
transfer calculations by changing the level populations. To 
quantify these excitation and ionization rates, the degrada- 
tion spectrum, which describes the number of fast electrons 
as a function of energy, needs to be known. There are several 
ways t o solve for the degradation spe c trum ( Spencer fc Fanol 
1954 : IShull fc van Steenberd Il985l: I Fransson fc Chevalierl 
19891 : IXu fc McCravl Il99ll : iKozma fc Fransson! Il992l ). One 



approach is the Monte Carlo method which, in general, is 
slow. Our approach is to solve the Spencer-Fano e quation 
following the method of iKozma fc Fransson! (|l992h . Once 
the degradation spectrum is known, all non-thermal excita- 
tion and ionization rates are computed and included in the 
statistical equilibrium equations (SEE). 

2.1 The Spencer-Fano Equation 

The radiative transfer code CMFGEN, which takes into ac- 
count radioactive decay, previously assumed that all radioac- 
tive decay energy is deposited as heat. In reality, fast elec- 
trons resulting from Compton scattering between 7-ray and 
thermal electrons lose energy through three different chan- 
nels - heating, excitation and ionization. The fractions of 
energy for the three channels can be computed if the degra- 
dation spectrum is known. This can be done b y solving 
the Spencer-Fano equat i on dSpencer fc Fanolll954T ). Follow- 
ing [Kozma~fcjiiransson| (|l992l ). the Spencer-Fano equation, 
which balances the number of electrons leaving and entering 
an energy interval, can be written as: 
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where A = min(_E max — E, E + I), iSmax is the maximum 
energy of non-thermal electrons, and I is the ionization po- 
tential. y(E)dE is defined as the flux of non-thermal elec- 
trons in the energy interval (E, E + dE). m is the number 
density of a single species, and the sum over i allows for 
the contribution of different ionization stages and species. 
L e (E) is the energy loss of the non-thermal electrons due 
to Coulomb scattering by the thermal electrons. Ej is the 
excitation energy, Oj is the excitation cross section of level j 
and a c (E,e) is the differential ionization cross section for a 
fast electron with energy E and an energy loss of e. There- 
fore, the ejected electron has energy e — I. The fast electron 
is called the primary electron and the ejected electron is 
called the secondary electron. The secondary electron can 
further degrade through the three same channels as the pri- 
mary electron. S(E) is the source term, and we assume all 
non-thermal electrons are injected at E max . The degrada- 
tion spectrum y(E) is normalized to £ max , so that y(E) is 
the degradation spectrum for 1 eV energy input. The first 
and second term on the left side represent the number of 
electrons leaving the energy interval (E, E + dE) by impact 
excitation and ionization, respectively. The term with L e (E) 
is the number of electrons leaving the energy interval by 
heating the medium. The first term on the right side is the 
number of electrons entering the energy interval (E, E + dE) 
by exciting another ion and losing energy. The second term 
is the ionization contribution from the primary electrons, 
and the third term is the ionization contribution from the 
secondary electrons if applicable. 

The Spencer- Fano equation is more conveniently solved 
nume rically using the integral form (|Kozma fc Franssonl 
EH), i.e., 

I>E r +E ' y(E')vj(E')dE' + y(E)L e (E) 

, _ /-Emax r{E'+E)/2 
+ Vn, / y(E') / a c (E' ,t)dtdE' 

Je Je'-e (2) 

= X>i / V(E') a c (E',e)dedE' 
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In Eq. [T] the integral of the differential ionization cross sec- 
tions needs to be done numerically. However, only the total 
ionization cross sections appear in the integral form of the 
Spencer- Fano equation (Eq. [2]), and can be computed ana- 
lytically. 

With the solution of the degradation spectrum, the frac- 
tional energies entering heating, rjh, excitation of ion i to 
level j, rjij, and ionization of ion i, r]i c , can be computed by 



^ = p— / y(E )L e (E )dE 
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where 



J2v(E + E^jiE + E l3 ) 
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+ / y(E')a tc (E',E + E)dE' 

J2E+Ii 

(Koz ma fc Franssonl Il992h . .Emit is the mean energy of the 
initial electrons, and Eo is the lowest excitation or ionization 
energy. 

To save computational effort, we only take into account 
the dominant ionization stages when computing the degra- 
dation spectrum. The other approximation is the neglect of 
7-ray transport. In reality, some 7-rays may escape the SN, 
while others may deposit their energy at large distances from 
where they were created. In these calculations we assume, 
for simplicity, in-situ deposition. A 7-ra y transport code 
has b een developed for use with CMFGEN (|Hillier fc Dessarq 
12012(1 

We previously mentioned that for the source term we 
assume that all high energy electrons are injected at E max . 
Although not entirely realistic, this kind of source function 
generally has little influence on the results. High energy elec- 
trons "forget" their initial energies after several scatterings, 
and the source function is quickly washed out. Alternatively 
we can assume a bell shape source function near E max . This 
tends to give a smoother degradation spectrum near E max . 



2.2 Ionization cross sections 



lArnaud fc Rothenflu j (|l985l ) adopted the formula proposed 
bv lYoungerl(|l98ll ) for calculating the impact ionization cross 
sections 



Qi(E) = -^{Ai 

ulf \ it 



Bi\l ) + 

u , 



+ Ciln(«) + A^} 



where u = E/E. E is the energy of the impact electron, 
and E is the ionization potential for the element. At, Bi, 
d, and Di are all coefficients. For the meaning of these 
coeffic ients, the reader should refer to lArnaud fc Rothenflua 
In our calculations, we have used the most up-to-date 
coefficients (Arnaud, private communication). 

With the total cross section computed by Eq. [7] the 
differential cross section a c (E,e) can be written as, 



a c (E, e) = Q(E)P(E, e — I) 



(8) 



E ,.< 



where P(E P ,E S ) is the distribution of t he secondary elec - 
trons for a primary electron energy of Ep. lQpal et all (|l97ll ) 
found P{E P , E s ) could be well described by 



© 2012 RAS, MNRAS 000,rflfT8l 



4 Chengdong Li, D. John Hillier & Luc Dessart 



P(E P ,E S ) = 



Jarctan[(£; p - 7)/2J] [1 + (£ S /J) 2 ] 



(0) 



where J is a parameter that gives a best fit to measurements 
and varies for different species. We adopt J = 0.6J, which i s 
shown to be a reasonable approximation (|Opal et al.lll97ll ). 
The secondary distribution function drops quickly so that 
very few secondary electrons have high energies. Specifically, 
for Hi with I — 13.6 eV, 66.7% of secondary electrons are 
below the ionization potential if E p = 1000 eV and 69.8% 
in the case of E v =10 000eV. Further, the mean energy of 
the secondary electrons, 7} a cc, only varies slowly as the en- 
ergy of the primary electron E v increases. For hydrogen, 
v scc = 21.53 eV for E v = 1000 eV and v scc = 33.36 eV for 
E p =10 000eV. Secondary electrons with low energies are 
easily thermalized and deposit energy as heat. 




4 6 

x=(mv 2 /2AE) 1 ' 



Figure 1. The tabulated gaunt factor (stars) for neutral atoms 
as a function of x ■ 



given in lvan Regemorte U dl962h . The 
second order (solid) and the third order (dashed) polynomial fit 



2AE 



arc also shown. 



2.3 Excitation cross sections 

The electron impact excitation cross sections are only avail- 
able for a few species, and only for a few levels. To pro- 
ce ed further, we adopt the Bethe approximation as discussed 
bv lvan Regemorterl (|l962l ). With this approximation, Qij is 
given by 



T, 



8tt 1 



Iff , _ 2 
Ji3 9 %a 



2 Ej-E 



(10) 



where fef is the initial electron kinetic energy normalized 
to hydrogen ionization energy 13.6 eV, Ih is the hydrogen 
ionization energy, Ej — Ei is the energy difference between 
the two levels, fij is the oscillator strength, and ao is the 
Bohr radi us. The Gaunt factor o for neutral atoms utilizes 
Table 1 of lvan Regemorte 3 (|l962h . which is shown in Fig. Q] 

The variable x is defined by x = yf^Jj where \ mv2 is the 
kinetic energy of the impact electron and AE is the transi- 
tion energy. For positive ions, g has a constant value of 0.2 
when x ^ 1 and is the same as the g for neutral atoms when 
x > 1. We also show the second order and the third order 
polynomial fit to the tabulated data. Although the third or- 
der gives a better fit to the data, it drops rapidly at larger 
values of x. Thus, we adopt the second order polynomial fit 
in our calculations. At high energies, the Bethe approxima- 
tion shows that the collisional excitation cross-section, like 
the ionization cross-section, scales as log E/E. 

Note that the Bethe approximation is only for permit- 
ted transitions (electric dipole transitions among states with 
the same spin). However, impact excitation by non-thermal 
electrons can occur for all transitions. To make allowance for 
those forbidden transitions, we use the collision strengths to 
compute the excitation cross sections, such that 



n 1 1 n ^ 



(ii) 



where Wi is the statistical weight of level i and is the 
collision strength. The main difficulty with this approach 
is that very few collision strengths are available. What we 
usually have are the thermally averaged effective collision 
strengths, defined by: 



-§V(§ 



(12) 



where Ej is the kinetic energy of the outgoing electron. For- 
tunately, the thermally averaged effective collision strengths 
for forbidden transitions are almost constant at high tem- 
peratures, and hence we adopt such constant effective colli- 
sion strengths as the collision strengths for impacting non- 
thermal electrons at all energies. This approximation is un- 
likely to have a major influence on the results since at high 
energies collisions for permitted transitions dominate (see 
Section I5.2|l , while at low energies there are relatively few 
non-thermal electrons, and collision rates are dominated 
by thermal electrons. An important distinction from non- 
thermal ionization is that non-thermal excitation does not 
produce secondary electrons. 

The accuracies and the uncertainties of the Arnaud & 
Rothenflug cross sections and the Bethe approximation will 
be discussed in Section \S. II 



3 THE HYDRO DYNAMICAL INPUT 

The hydrodynamical model l lml8a7Ad' (Woosley, private 
communica tion) we use as an i nitial input is the same as that 
adopted bv lDessart fc Hillierl |2010l ). We briefly summarize 
the main properties of this model below. The hydrodynami- 
cal m odel was produced by the code kepler (| Weaver et al.l 
Il978h . using a progenitor with main-sequence mass of 18 
M© . The metallicity of the Large Magellanic Cloud (LMC) 
was adopted, which is Z = OAZq. The star is a BSG when 
the explosion is initiated. 

The hydrodynamical model at 0.27 day was employed, 
since at this time homology is a good approximation in the 
outer regions. The rest of the ejecta is enforced to be strictly 
homologous to accommodate CMFGEN, although the inner- 
most part of the ejecta is trimmed since it suffers fallback. 
Moderate mixing is induced manually to soften the compo- 
sition gradients. Hydrogen is deficient below 1500 kms -1 , 
while 56 Ni is primarily found below 2000 kms -1 . The to- 
tal amount of 56 Ni in the initial hydrodynamical model is 
0.084 Mq. Although a large set of model atoms are utilized, 
some potentially important species such as Til, Fel, Scl 
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Figure 2. Illustration of the elemental mass fractions on a log- 
arithmic scale as a function of velocity at day 127. We plot hy- 
drogen (black), helium (red), carbon (green), oxygen (blue) and 
56 Co (magenta) — five important elements in the model. 

and Sell are missing. The pre-SN steady and explosive nu- 
cleosynthesis in the kepler model were undertaken with a 
13-isotope network which only captures the composition ap- 
proximately. It also only computes explicitly the abundance 
of one unstable isotope ( 56 Ni) - numerous intermediate-mass 
elements (IMEs) and iron-group elements (IGEs) are bun- 
dled together. Moreover, the hydrodynamical model is en- 
forced to be homologous, which affects the kinematics of the 
inner ejecta. With the recent availability of a model atom 
for Fel we discuss the influence of Fel on the spectra in 
Section [()] 

The models from day 0.27 to day 20.8 have been pre- 
sented in iDessart fc HillieU (|201Cf ) . focusing on early evolu- 
tion at the photospheric phase. 



4 THE NON-THERMAL MODEL 

To study the long-term evolution of our SN model, we 
evolved the model further in time until day ~ 1000 - well 
into the nebular phase. After about 40.6 days, the sequence 
is separated into two branches, a "thermal sequence" with 
7-ray deposition as heat exclusively, and a "non-thermal se- 
quence" with non-thermal excitation and ionization. 

In this section, we focus on a typical non-thermal model 
and discuss the properties of non-thermal processes in the 
model. We select a non-thermal model at day 127 (here- 
after model D127_NT) near the beginning of the nebular 
phase. It is during the nebular phase that non-thermal ef- 
fects will be the most important. Moreover, complications 
arising from dust will not be present s ince dust formation 
is not thought to happen until day 300 (|Suntzeff fc BouchetJ 
Il99d : ICehrz fc Nevlll990f ) . 

In Fig. [2] we illustrate the elemental mass fractions on 
a logarithmic scale as a function of velocity at day 127. Hy- 
drogen is an order of magnitude less abundant at velocity 
1500 kms -1 and this velocity can be taken to indicate the 
velocity down to which H has been mixed. Rather than 56 Ni, 
we show the mass fraction of 56 Co, since at day 127 only a 
small amount of 56 Ni is left, and 56 Co decay is the main en- 
ergy source. The figure shows that 56 Co has been mixed out 
to about 2000 kms" 1 in the initial model, although trace 
amounts also exist at higher velocities. 
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Figure 3. The degradation spectra at a model velocity of ~ 
1000 kms -1 computed using different numbers of energy bins: 
100 (black), 1000 (red), 10 000 (blue). 



4.1 The degradation spectrum 

A typical degradation spectrum at velocity ~ 1000 kms -1 
for the D127_NT model is shown in Fig. [3] (red). The sharp 
rise at high energy is a result of the adopted source function 
- we inject all fast electrons at E max . The spectrum from 
E max = 1 keV down to 100 eV shows the slowing down of 
primary electrons, while the rise at low energy indicates a 
cumulation of secondary electrons. The shape of the degra- 
dation spectrum, particularly at low energies, is influenced 
by the ionization state of the gas. 

Uncertainties in the degradation spectrum are intro- 
duced by the choice in the number of energy bins, accu- 
racies of the excitation and ionization cross sections, and 
by the choice of the high energy cut, E max . The latter two 
will be discussed in Section [9] For the number of energy 
bins, we typically adopt N = 1000 extending from leV to 
1 keV. In Fig. [3] we also show another two degradation spec- 
tra computed with energy bins N = 100 (black) and N = 
10 000 (blue). The flux in the N = 100 degradation spec- 
trum is systematically lower than the other two, except at 
the high energy edge. The N = 10 000 spectrum is very close 
to the one with N = 1000, and hence the uncertainty aris- 
ing from using 1000 energy bins is negligible. We also made 
spectral comparisons among the above three models. The 
models with 1000 and 10 000 energy bins showed no distin- 
guishable differences, while the model with 100 energy bins 
was only slightly different from the other two. Interestingly, 
we later discovered that a smaller number of bins can be 
used, provided we scale the degradation spectrum so that 
the fractional energies entering the three channels sum to 
unity. 

4.2 Number density of non-thermal electrons 

Knowing the degradation spectrum y(E), the electron spec- 
tral density in space can be estimated from y(E)/v, where 
v = \/2E/m e . The number density of the non-thermal elec- 
trons is then given by 

N ,=l--m dE =i'--„ im ^ dE (I3) 

Fig. [4] shows the comparison of number density of the 
non-thermal and thermal electrons at day 127. The num- 
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Figure 4. Comparison of the number density of the non-thermal 
electron (solid) and the thermal electron (dashed). 
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Figure 5. The fractional energy that goes into heating 5i?h(red), 
ionization SEi (blue) and excitation 5E e (green) in model 
D127_NT, as a function of velocity. The ejecta ionization fraction 
X e , which is defined as the ratio of number of ions to number of 
atoms, is also shown (black). 



ber density of the non-thermal electrons is many orders of 
magnitude smaller than that of the thermal electrons at all 
depths. The increase of the non-thermal electrons number 
density at inner regions is consistent with the increase of the 
Co abundance shown in Fig.fJ 



4.3 Energy fraction of the three channels 

With the inclusion of non-thermal excitation and ioniza- 
tion, 7-rays deposit energy into all three channels. In Fig. 
we plot the fraction of energy entering the three channels 
as a function of velocity. We also plot the ejecta ioniza- 
tion fraction X e , which shows that the heating fraction is 
sensitive to X e . In the outermost region, the heating frac- 
tion reaches unity as X e becomes 1 - the highest ionization 
stage ions are unable to be ionized and excited, quenching 
the excitation and ionization processes. The high ionization 
fraction in this region is set by ionization freeze-ou t due 
to the time-dependent effects (|Dessart fe Hi llier 2008). Be- 
tween 3000 kms" 1 and 10 000 kms" 1 , non-thermal excita- 
tion and ionization become prominent where X e becomes 
very small. As X e increases inward from 3000 kms -1 , frac- 
tional non-thermal excitation and ionization decline and 
fractional heating rises and becomes dominant. 



Figure 6. The fractional energy (solid lines) consumed by a 
species in the non-thermal excitation (top) and ionization (bot- 
tom) processes as a function of ejecta velocity. Only the largest 
5 (summation over all depths) fractions are shown. The fractions 
of species abundances are plotted in dashed lines for comparison. 
Importantly, the highest ionization stage of each element is ex- 
cluded in the calculations of fractional abundances, (see the text 
for details). 



4.4 Excitation and ionization 

When non-thermal excitation and ionization occur, energy 
from non-thermal electrons is used to excite or ionize an 
atom. The Spencer-Fano equation (Eq. [1} implies that the 
abundances and the cross sections are the key to determine 
what energy a species "consumes". In this subsection, we 
study the fractional energies consumed by a species and 
compare them to the corresponding species. 

In Fig. [6] we illustrate the fractional energies taken away 
from the non-thermal electrons by a species in both the ex- 
citation and ionization channels. The largest 5 consumers, 
summing over all depths, are plotted (solid lines). We over- 
plot the fractional abundances of these 5 species (dashed 
lines). However, the calculations of fractional abundances 
excludes the highest ionization stages of all elements. These 
highest ionization stages are the ones in the model, not nec- 
essarily the same as the ion of an element with no electron 
left. For instance, in our models, we utilize detailed model 
atoms for Fe II, Fe in, Fe iv, Fe v, Fe vi and Fe vn, then Fe vm 
is the highest ionization stage for iron. The reason for ex- 
cluding the highest ionization stages is that we treat these 
ionization stages as a single level, and so they don't have 
any excitation or ionization routes. In both channels, it is 
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Figure 7. Comparison of temperature structure between model 
D127.NT (solid) and model D127 (dashed). 
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obvious that Hi, He I and Hell are the dominant consumers 
of non-thermal electron energy. Hell is the main consumer 
above 30 000 km s -1 . In the middle region, both H I and He I 
are prominent. Due to low abundance of H, He I dominates 
the inner zone, although other species, like Cl, Cn, Ol and 
O II also play a role. Fig. [6] shows that the species abun- 
dances, especially Hi, He I and Hell, have similar patterns 
to their fractional energies, in both channels. Such similarity 
indicates that allocation of the non-thermal energy is gen- 
erally controlled, as expected, by the species abundances. 
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5 NON-THERMAL VS THERMAL MODELS 

In this section, we present a comparison of the temperature 
and ionization structure between the non-thermal and ther- 
mal models. The spectral evolution of the two sequences will 
also be compared. 

5.1 Comparison of the temperature structure 

In Fig. [7J we compare the temperature structure between 
model D127.NT and the thermal model at day 127 (here- 
after model D127). As mentioned previously, the non- 
thermal model generally has a lower temperature in the re- 
gion between 1000 kms -1 and 5000 kms -1 where the non- 
thermal effects are crucial. The temperature for the two 
models above 5000 kms -1 is pretty flat (i.e., almost con- 
stant). This is an artifact - a floor temperature of 1500 K is 
imposed for the models to prevent numerical overflow. The 
flux mean opacity at 5000 kms -1 is of order 10 -2 , thus the 
region above 5000 kms -1 has very little influence on spectral 
formatiorQ. 

The decrease in the temperature can be understood as 
follows. The thermal model puts all decay energy into heat- 
ing, while the non-thermal model only puts a portion of 
decay energy into heating. Moreover, the coupling between 
the gas and the radiation at this time is relatively weak, 
and hence, we expect that only the heating channel will sub- 
stantially affect the temperature. Energy that has not gone 
into heat/temperature is stored up in ionization/excitation, 

1 This is confirmed by computing two models with different floor 
temperatures. The synthetic spectra show no distinguishable dif- 
ference between models with floor temperatures of 1500 K and 
4000 K. 



Figure 8. Logarithmic species fractions of different ionization 
stages for hydrogen (top) and helium (bottom). Comparison is 
made between model D127_NT (solid lines) and D127 (dashed 
lines). 

driving the material even further away from LTE level pop- 
ulations. 

Notice that the temperature of the two models is simi- 
lar below 1000 kms -1 , despite the existence of a substantial 
amount of 56 Ni. Even though the region is not in LTE, colli- 
sional processes are important, and there is a strong coupling 
between the radiation field and the gas, and thus energy ini- 
tially deposited via the ionization and excitation channels 
can be thermalized. The fractions of energy entering the 
three channels are shown in Fig. [5] The heating fraction be- 
comes dominant in the inner region, which is consistent with 
the temperature comparison in Fig. [7] 

5.2 Controlling processes for Hi and Hel lines 

In this section, we compare the hydrogen and helium ioniza- 
tion structure in the D127 and D127_NT models and study 
the processes controlling the Hi and Hel lines in the non- 
thermal model. The H and He ionization fractions for both 
the D127 and D127.NT models are illustrated in Fig. [8] In 
the cobalt region, the ionized H fraction has increased by up 
to an order of magnitude, while in the region above v > 5000 
kms -1 , the ionization is unchanged, since 56 Co is deficient 
there. 

While He remains neutral for v < 10 000 kms -1 , the 
singly ionized He fraction is now significant. We can also 
see significant change in the He ionization fraction above 
the region where cobalt is abundant, which is due to a 
small amount of cobalt in the model. At large velocities, 
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Figure 9. Upper panel: The D127 model spectrum (black) and 
the model spectra computed by including only bound-bound tran- 
sitions of He I (red) and Hi (blue) at optical band. Lower panel: 
Same as the upper panel, but for the D127-NT model. 
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Figure 10. Illustration of model spectrum (solid) and model 
spectrum computed by including only bound-bound transitions 
of He I (dashed) for model D127 (top panel) and D127.NT (bot- 
tom panel) between 1.9 and 2.2 /im. 



the ionization of both H and He is frozen, since the expan- 
sion time-scale is smalle r than the recombination time-scale 
jDessart fc HillierlbOldl ). 

Fig. [9] shows the model spectra and the spectra with 
only bound-bound transitions of Hel and Hi at the opti- 
cal band. Optical, as well as IR, Hi and Hel lines mainly 
originate below 2000 kms -1 in the non-thermal model, 
with a small fraction of the line flux coming from 2000 - 
3000 kms -1 for some lines, e.g., Hel 1.083 ^m. Hi lines in 
the thermal model (upper panel) are generally weak and He I 
lines are absent. In the non-thermal model (lower panel), all 
Hi lines are strengthened, with Ha particularly boosted by 
non-thermal processes. Hel lines also contribute significantly 
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Figure 11. Top: Normalized rates for different processes that 
populate and depopulate the ls2p 1 P° state of Hel. Only pro- 
cesses with a fractional rate ^ 5% at more than one depth 
are shown. Bottom: Same as the top panel, but for the case of 
ls2p 3 P° state of Hel. 



to the spectrum. While most optical Hel lines are contami- 
nated by other lines, Hel 7065 A might be a good diagnostic 
to infer the influence of non-thermal processes on helium. 
Notice that He I 4471 A is excited in the non-thermal model 
but is absent from the spectra due to line- blanketing. 

The prominence of Hel 1.083 (im and Hel2.058/mi 
due to non-ther mal processes at lat e times (t ^ 200 
d) was noted by iLi fc McCravl (| 19951), and observation s 
showed He 1 2.058 /im is more isolated (jMeikle et aHll989h . 
In Fig. I lOt we show the model spectra and the spectra with 
only bound-bound transitions of Hel for models D127_NT 
and D127, focusing on the He 1 2.058 /im line. Similar to 
the optical band, no Hel line is present in model D127 
in the wavelength range from 1.9 /im to 2.2 /im. However, 
He 1 2.058 /im produces a prominent absorption feature in 
the D127_NT model. This feature has b een seen in obser- 
vations of SN 1987A (jMeikle et al.lll989T ). providing strong 
evidence for t h e influ ence of non-thermal processes on He I. 
ILi fc McCravl (|l995h showed that He 1 2.058 (im is mainly 
due to non-thermal excitation and He 1 1.083 /mi is mainly 
excited by thermal electrons during 200 d ^ t ^ 450 d. 

In Fig. 1111 we illustrate the major fractional rates 
5% at any depth) that populate and depopulate the states 
ls2p 1 P° and ls2p 3 P° of Hel for model D127.NT. Posi- 
tive normalized rates illustrate the routes that populate the 
state, while negative ones show the routes that depopulate 
the state. States Is 2p 1 P° and Is 2p 3 P° are the upper levels 
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Figure 12. Normalized rates for Hi (top) and He I (bottom), 
including all processes allowed in the modeling. 



for the He I lines 2.058 /im and 1.083 /im, respectively. For 
the state ls2p 1 P°, we confirm that non-thermal excitation 
is the main excitation mechanism. Cascades from higher lev- 
els, especially the ls3d 1 D state, also play an important role 
at this stage. 

For the state ls2p 3 P°, cascades from higher lev- 
els (mainly ls3s 3 S and ls3d 3 D) are the key processes 
to populate it below 3000 kms -1 . Recombination is also 
non-negligible. These are indirectly related to the non- 
thermal ionization processes. Thus, in contrast to the 
He 1 2.058 /^m line, non-thermal excitation is largely irrel- 
evant for He 1 1.083 /^m. This is due to the difference in 
the collision strength behavior of the two transitions. Non- 
thermal excitation from the ground state (Is Is S) is gen- 
erally the dominant excitation route, since the ground state 
is the most populated state. At low energies, the collision 
strength for the singlet transition (permitted transition), 
SI (Is Is L S - Is 2p 1 P°), is similar to that of the triplet transi- 
tion (forbidden transition), Q(ls Is X S - Is 2p 3 P°). However, 
f2(ls Is X S - Is 2p 1 P°) grows significantly as the impact en- 
ergy increases, while $7(ls Is X S - Is 2p 3 P°) is approximately 
constant at high energies. As a result, non-thermal excita- 
tion processes contribute much less to the triplet states than 
to the singlet states. Due to radiative excitation from Is 2s 3 S 
to ls3p 3 P° and subsequent decay to ls3s 3 S, electron cas- 
cade from ls3s 3 S, rather than ls3d 3 D, is the dominant 
process populating the ls2p 3 P° state above 3000 kms -1 . 

Despite the importance of non-thermal processes, we 
find that photoionization and recombination play a crucial, 
if not dominant, role in populating many other levels. How- 
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Figure 13. Comparison of optical and IR spectra between models 
D127.NT and D127. 
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Figure 14. Comparison of the Ho profile between models 
D127.NT (solid) and D127 (dashed). 



ever, in a "thermal" model, the photoionization and recom- 
bination rates are too small to produce a prominent effect. 
The trigger for the large photoionization and recombination 
rates is non-thermal processes. In Fig. [12] we show the nor- 
malized rates affecting the H I and He I ionization balance for 
model D127_NT. For Hi, the fraction of non-thermal energy 
deposited via the ionization channel is small and photoion- 
ization and recombination are the dominant processes and 
they balance each other below 8000 kms -1 . We check an 
individual level, n = 1, of Hi in both models (D127_NT and 
D127), and find that the photoionization and recombination 
rates are increased by a factor of ~ 2 when non-thermal pro- 
cesses are included. For the He I Is 2p 1 P° state, non-thermal 
ionization and excitation rates in model D127_NT are 5 or- 
ders of magnitude greater than any rate in model D127. 
Below 3000 kms -1 , non-thermal ionization becomes signifi- 
cant, and the radiative recombination rate is now balanced 
by the photoionization and the non-thermal ionization rates. 
Although the net non-thermal rate is not overwhelmingly 
dominant, it is crucial for controlling the ionization balance. 
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Figure 15. Left panel: montage of synthetic optical spectra for the thermal sequence. The days since breakout are labeled below each 
spectra. Right panel: montage of synthetic spectra at the same epochs shown in the left panel for the non-thermal sequence. The synthetic 
spectra are reddened with E(B-V) = 0.15 and scaled for a distance of 50 kpc. The most prominent difference between the spectra of the 
two sequences is the evolution of Ha - its strength persists at all times in the non-thermal sequence but it almost disappears at late 
times in the thermal sequence. 



5.3 Comparison of optical and IR spectra 

Fig. [13] shows the comparison of the optical and IR spectra 
between models D127_NT and D127. In the optical range 
(Fig. 1131 top panel), D127_NT generally produces slightly 
lower fluxes but with two enhanced prominent features — 
Ha and He I 7065A. The Can IR triplet is slightly weakened. 
In the bottom panel of Fig. 1131 most of the lines are strength- 
ened, such as He 1 1.083 /im, Pa 7 1.094 fim, Ol 1.129 /im, Pa 
/3 1.281 fim and Pa a 1.875 fim. However, Call 1.184 fim and 
1.19 5 ^m are weakene d. 

iDessart fc Hillier! (|201lh found that in their SNe II-P 
models the Balmer lines suddenly disappeared at end of the 
photospheric phase - a result of the vanishing of Balmer- 
continuum photons. However, observations of Type II SNe 
show a strong Ha profile during the nebular phase. They 
attributed this discrepancy to the exclusion of non-thermal 
excitation and ionization in the models. We present here the 
comparison of Ha at day 127 in Fig. 1141 While Ha is almost 
absent in model D127, it shows a strong P Cygni Ha profile 
in model D127_NT. The emission component is stronger and 
wider and a P Cygni absorption trough is now present. The 
width of Ha is related to the outward mixing of 56 Ni. This 
may be helpful to determine the 56 Ni mixing in Type II SNe, 
as hydrogen is abundant in these SNe and Ha is one of the 
strongest features at nebular epochs. 



5.4 Comparison of the spectral evolution 

As m entioned in the previous section, IDessart fc Hillierl 
|201ll ) were unable to reproduce the strong Ha at the nebu- 
lar phase for red supergiant (RSG) models, which was asso- 
ciated with the neglect of non-thermal processes. Moreover, 
the mixing in their models may have been too weak. The 
best explanation for the persistence of Ha at late times is 
non-thermal excitation and ionization processes combined 
with mixing of H and 56 Ni. This has been one of the mo- 
tivations for this study. In Fig. 1151 we present the spectral 
evolution of the non-thermal and the thermal sequences. The 
disappearance of Ha is also observed in the thermal model 
for SNe resulting from a blue supergiant (BSC) progeni- 
tor. The strength of Ha increases till day 96, and then fades 
quickly. On day 140, no prominent feature is seen. The spec- 
tral evolution of the non-thermal sequence shows similarity, 
except for Ha, which has a stronger maximum strength at 
about day 96, decreases gradually after that, and maintains 
a stronger signature on day 140. This confirms that the per- 
sistence of strong Ha emission during the nebular phase is 
closely related to non-thermal excitation and ionization. The 
non-thermal effects on spectra at day 54 are relatively small, 
confirming that the neglect of the non-thermal processes in 
earlier models is unlikely to have a significant influence on 
model spectra. 
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Figure 16. Comparison of optical spectral between 
model D127.NT and the model with Fel (hereafter model 
D127.NT.FcI). 
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Figure 17. Comparis on of the observed (black stars, 
ISuntzeff fe Bouch ct 1990) bolometric light curve of SN 1987A to 
the theoretical bolometric light curves of the thermal (red) and 
non-thermal (blue) sequences, with symbols referring to the com- 
puted epochs of models. The observational bolometric light curve 
is constructed by using ultraviolet, optical and infrared photome- 
try. A light curve resulting from the conversion of energy from the 
radioactive decay of 0.084 Mq (green) is also shown for reference. 



6 THE INFLUENCE OF Fe I 

The source of opacity in SNe is complicated, since lots of 
IMEs and IGEs are synthesized during the explosion. Al- 
though a huge number of levels are included in the model, 
we are still missing some model atoms and the amount 
of levels of current model atoms may also be insufficient. 
Many of the missing model ions are the lowest ionization 
stages. These do not influence the spectra at early times 
when the ejecta is hot, but may significantly influence the 
opacity at late times when the ejecta becomes cold. Here, 
we e xplore the in fluence of Fel. Oscillator strengths are 
fromlKuru cz ( 2009^3: photoionization cross sections are from 
iBautistal (| 19971 ) (accesed through Nahar-OSU-Radiative- 
Atomic-Data web site), while collision rates among the low- 



2 Data from R. L. Kurucz 
http://kurucz.harvard.edu 
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Figure 18. Comparison of the V-, R-, and I-band light curves 
of SN 1987A to the corresponding synthetic light curves. The 
observed V, R, and I magnitudes are shown in stars, and the 
synthetic V, R, and I magnitudes from the non-thermal sequence 
arc plotted in squares connected by dashed-dottcd lines. A scaling 
is applied to the R and I band magnitudes, both observed and 
synthetic, to optimize visualization. The synthetic photometry is 
computed by con volving the sy nthetic spectra with the Landolt 
filter bandpasscs iLandohj|l9"92[) . We applied a reddening of E(B- 
V)=0.15 and adopted a distance of 50 kpc. 



est 10 levels are from lPelan fe Berringtonl (|l997l ) (accessed 
from T1PBASE at cdsweb.u-strasbg.fr/tipbase/home.html). 

We included Fe I at day 5 and ran a sequence of models 
with this new atom and with non-thermal excitation and 
ionization. Fig. [16] illustrates the comparison of the optical 
spectra at day 127. With the inclusion of Fel, the emergent 
flux decreases appreciably between 3000 A and 5500 A, with 
the energy redistributed to longer wavelengths above 6000 
A. Many other atoms, such as Col, Til and Sell, may cause 
effects similar to Fel. We will include more complete model 
atoms, as they become available, to better address line blan- 
keting. 



7 COMPARISON WITH THE OBSERVATIONS 

SN 1987A is still one of the best observed and studied 
SNe, and it is about 50 kpc away from the earth. Al- 
though its progenitor was a BSC, there is still debate 
about the evol utionary channel that l ed to a SN explo- 
sion as a BSC (iHillebrandt et al.1ll987l ;l Langer et al. Hl989l ; 
IPodsiadlowski et al.lll99ol ). In this section, we present photo- 
metric and spectroscopic comparison with the observations 
of SN 1987A. Our observational data is taken from CTIO 
(jPhillips et al .11988^ For the synthetic sp ectra, we adopt the 

and use a redden- 



extinction curve of 



Cardelli et al 



LCtlC spi 

il989|) 



ing of E(B-V)= 0.15. iPanagial (|2005l ) derived the distance 



to SN 1987A to be 51.4 kpc with an unce rtainty of 1.2 kpc. 
We ad opt d = 50 kpc for consistency with iDessart fe Hillierl 
(|2010t ). 
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Figure 19. Comparison between observed spectrum of SN 1987A and the model spectra at day 127. A scaling is applied to the observed 
spectrum to correct for the difference between the spectroscopic (4.30) and the photometric (4.37) V-magnitude. The three model spectra 
are scaled by 0.07/0.084 according to the 56 Ni mass in SN 1987A and in the model. 



7.1 The synthetic and observed light curves 

The comparison of the bolometric light curves is illustrated 
in Fig. 1 171 while Fig. 1 181 illustrates the comparison of the V-, 
R-, and I-band light curves. The theoretical bolometric light 
curve rebrightens at a later time than observation, which is 
probably due to insufficien t outward mixing of 56 Ni |Utrobinl 
120041 : iDessart et al.ll20l3 ). The delay of the predicted peak 
luminosity supports this idea. Both the thermal and non- 
thermal bolometric light curves predict a higher luminosity 
than th e observed one, partly due t o higher 56 Ni mass in our 
models. ISuntzeff fc Bouchetl (|l990l ) derived a 56 Ni mass of 
0.07 M Q , while we have 0.084 M Q in the model, which is the 
amount that comes out of the explosion model. As the pur- 
pose of this work is to explore the influence of non-thermal 
processes, and not to fit the spectra of SN 1987A, we did not 
adjust the progenitor model, nor did we explore the influence 
of mixing, explosion energy, and mass cut in the progenitor 
model. The difference in the amount of 56 Ni only affects 
our results quantitatively, not qualitatively. The hydrody- 
namical model, such as the size of the helium core and the 
enforcement of homology at the beginning of the modeling, 
also introduces uncertainties. Considering that we have no 
free parameter, the synthetic bolometric light curves agrees 
reasonably well with the observed bolometric light curves. 
We also plot a light curve resulting from the radioactive de- 
cay of 0.084 M of 56 Ni. The agreement between this light 
curve and the synthetic bolometric light curve during the 
nebular epoch demonstrates that the energy source of the 
ejecta is totally dominated by radioactive decay. 

Although the synthetic bolometric light curve is reason- 
ably well reproduced, the synthetic multi-band light curves 
show varying degrees of systematic offsets. We are mainly in- 
terested in the nebular epochs when the impact of mixing is 
less important. The R-band magnitudes are close to the ob- 
servations ~ 130 days after explosion. However, the V- and 
I-band magnitudes show significant discrepancies (Fig. I18[) 
on the tail of the light curve, with the V-band overestimated 
by ~ 1 magnitude and the I-band underestimated by ~ 0.4 
magnitude. The difference in the 56 Ni mass between the 
models and SN 1987A translates into ~ 0.2 magnitude in 
all bands. This does not remove the discrepancy in the V 
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Figure 20. Same as Fig. 1191 but it now shows the comparison of 
the Ha profile in the velocity space between observation (solid) 
and model D127.NT (dashed). 



and I bands. However, the overestimation in the V-band 
and the underestimation in I-band indicate that strong line- 
blanketing effects are missing in the models. These would 
redistribute the V-band fluxes to the I-band, possibly alle- 
viating the discrepancy. Future studies with more complete 
model atoms may quantify such redistribution. In fact, the 
non-thermal sequence with Fel reduces the V magnitudes 
by ~ 0.2 magnitude. 

7.2 Spectral comparison 

Fig. [T5] shows a comparison of the observed optical spec- 
tra against the model spectra of D127.NT (red), D127 
(blue) and D127_NT_FeI (green). The Ha profile is too weak 
in model D127, while both non-thermal models, D127.NT 
and D127_NT_FeI, show considerably stronger Ha. Mod- 
els D127_NT and D127 show much stronger fluxes between 
3000 A and 6000 A, which is reflected by the 1-magnitude off- 
set in the V-band magnitudes in the previous section. The 
inclusion of Fel improves the fit, mainly in the wavelength 
range from 3000 A to 4500 A. 

In Fig. 1201 we show the theoretical and observed Ha 
profiles in velocity space. While the peak of Ha in the non- 
thermal model is comparable to that of the observation the 
profile is narrower. The left emission wing of the profile fits 
the observation, but the Ha profile in the observation ex- 
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tends at least 3000 kms -1 further to the red than does the 
theoretical profile. Further, the P Cygni absorption compo- 
nent does not match — the minimum of the model absorp- 
tion sits at 3000 kms -1 , while observational minimum is at 
5000 kms" 1 . 

The full width at half maximum (FWHM) of Ha in the 
synthetic spectra is about 4000 kms -1 , indicating that the 
main region that contributes to Ha is the core zone below 
2000 km s - 1 . This is coincident with the region in which 56 Ni 
is abundant. The spatial distribution of 56 Ni depends on the 
degree of mixing. Mixing effects in SNe explosion are widely 
seen in 2D and 3D models, however in pure ID hydrody- 
namic models mixing must be induced artificially. According 
to the observed broad Ha profile, our ID hydrodynamical 
model has insufficient mixing of 56 Ni. To constrain the mix- 
ing is not the goal of this paper, but the influence of mixing 
will be discussed in Section [9] Another possible explanation 
is 7-ray transport. We made the assumption that all 7-rays 
deposit energy locally, but 7-ray transport makes it possible 
to deposit more energy at larger velocities. A discussion of 
the possible influence of 7-ray transport will be carried out 
in Section [9] 



7.3 Comparison of the spectral evolution 

To better compare the spectral features between the obser- 
vations and the synthetic spectra, we selected the synthetic 
spectra which are about 10 days later than the observations. 
In Fig. 1211 the observations are shown on the left panel and 
the synthetic spectra are shown on the right panel. Recall 
that Balmer lines disappear suddenly at the end of the pho- 
tospheric phase in the thermal sequence, thus the biggest 
improvement of the non-thermal sequence is the persistence 
of the strong Ha profile at late times. In the non-thermal 
sequence, the strength of the Ha profile first increases, and 
then decreases slowly, which is consistent with the obser- 
vations. The NalD lines are contaminated by He 1 5875 A. 
However, NalD lines show behavior similar to Ha, being 
narrower and weaker than the observations. This is also re- 
lated to the non-thermal effects, since higher H ionization 
produce s more electrons, which leads to an increase in neu- 
tral Na jDessart fc Hillierll2008h . The synthetic spectra are 
rich with weak lines around 4000 A, while the observations 
show few weak lines there. Fen is the main contributor to 
these lines. With the inclusion of Fel, this discrepancy be- 
comes much smaller. Moreover, the model does not contain 
scandium and barium, which may explain a few missing 
P Cygni profiles in the synthetic spectra. The absence of 
Ban 6142 A is suggestive. The absence of Sell in our model 
likely le ads to an underestimat e of line blanketing around 
4000 A (|Dessart fc Hillieill201lf) . 



8 UNCERTAINTIES 

8.1 Impact excitation and ionization cross 
sections 

The uncertainty of the non-thermal solver mainly comes 
from the uncertainties of the electron impact excitation 
and ionization cross sections. As mentioned previously, the 
impact excitation cross sections are computed with the 



Table 1. Combinations of species to scale the cross sections 



Model 


excitation 


ionization 


EH5 


H and He by 5 


Same 


EM5 


Metal species by 5 


Same 


EA5 


All species by 5 


Same 


IH5 


Same 


H and He by 5 


IM5 


Same 


Metal species by 5 


IA5 


Same 


All species by 5 



Bethe approximation, and the impa ct ionization cross sec- 
tions utilize the analytical formula of lArnaud fc Rothenflud 
|l985l ) with updated coefficients. The Bethe approximation 
is a high-energy approximation derived from the Born ap- 
proximation that neglects the short range interaction be- 
tween the perturbing electron and the atomic electron. 
Ivan Regemorterl (|l962 > corrected the Gaunt factor accord- 
ing to existing experimental data and accurate calcula- 
tions to allow the formula to be used at low energies. The 
Bethe approximation generally works better at high en- 
ergies. For some less understood species, the Bethe ap- 
proximation may brea k down even at high energies. The 
lArnaud fc Rothenflud formula may have troubles for ele- 
ments heavier than helium. Both experimental and theoreti- 
cal cross sections are required to evaluate their u ncertainties, 
but both of them are very difficult to estimate |Burke et al.l 
1 1994 iRamsbottom et alil2007l ; lRamsbottomll2009l '). 

To investigate the influence of the uncertainties, we 
ran six test models by artificially scaling the impact exci- 
tation or ionization cross sections for different combinations 
of species. These test models were all based on the same 
input model computed using standard cross-sections, and 
were run until the populations and the radiation field were 
fully converged. Table [1] lists the combinations and scalings 
applied to the cross sections. In Fig. 1221 we compare the 
fractional energies of the three channels in these six tests. 
In the case of varying the excitation cross sections (Fig. 1221 
left panel) , increasing the excitation cross sections of H and 
He results in an increase of fractional energy in the excita- 
tion channel and a decrease of fractional energy in the ioniza- 
tion channel. The change in the heating fraction is relatively 
small. An increase in the excitation cross sections of metal 
species has no effect above 2000 kms - , which is obviously 
due to their abundance deficiency. Below 2000 kms -1 , this 
also tends to increase the fraction of energy in the excitation 
channel and to decrease the fraction of energy in the ioniza- 
tion channel. The ejecta ionization fraction X e in all cases 
is roughly the same, implying that the ionization structure 
of the ejecta is hardly affected. Varying the ionization cross 
sections (Fig. 1221 right panel) has the reverse effect of vary- 
ing the excitation cross sections - the increase of ionization 
cross sections results in an increase in the fractional energy 
of the ionization channel and a decrease in the fractional en- 
ergy of the excitation energy. This is what we expect from 
the Spencer- Fano equation, given that the source term S(E) 
is fixed. 

The fractional heating remains almost unchanged be- 
tween 2500 kms -1 and 10 000 kms -1 when the ionization 
cross sections increase. This surprising behavior is related to 



a low ejecta ionization fraction, X e , seen in Fig. 22(b) At 
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Figure 21. Left panel: montage of observed optical spectra of SN 1987A. The days since breakout are shown below each spectra. Right 
panel: montage of synthetic spectra of the non-thermal sequence. Each observed spectrum on the left is compared to a synthetic spectrum 
at roughly the same time since explosion. The synthetic spectra are reddened with E(B-V) = 0.15 and scaled for a distance of 50 kpc. 
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Figure 22. Comparison of the energy fractions of the three 7-ray decay channels in models in which we have artificially scaled the 
excitation or ionization cross sections. The model used for testing is the one at day 127 without including Fel. Left: The 3 test cases of 
scaling the impact excitation cross sections. The scaling conditions are denoted in the figure. Different colors arc fractions of different 
channels. The solid lines shows the fractions of the original model. The dotted lines, the dashed lines and the dashed-dotted lines represent 
the EH5 model, the EM5 model and the EA5 model, respectively. Right: The same as the left, but it shows the 3 cases of scaling the 
impact ionization cross sections. The dotted lines, the dashed lines and the dashed-dotted lines show the channel fractions of the IH5, 
IM5 and IA5 model, respectively. 
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these low-ionization fractions, the degradation spectra show 
that almost all electrons have degraded into very low ener- 
gies and can no longer cause ionizations. Since an enhance- 
ment in the ionization cross-section does not dramatically 
change the number of secondary electrons produced, the ion- 
ization channel gains energy at the expense of the excitation 
channel. Conversely, a change in the excitation cross-sections 
does change the amount of energy entering the thermal chan- 
nel. Crudely, excitation causes higher energy losses before 
the creation of a secondary electron via ionization. The com- 
petition between various processes is extremely sensitive to 
X e . As X e becomes appreciable (e.g., above 10000 kms -1 .), 
many degraded electrons possess high energy. Increasing ei- 
ther the ionization or the excitation cross sections strength- 
ens non-thermal processes at the expense of thermal energy 
deposition. 

There are two types of uncertainties introduced by the 
approximations in computing the excitation and ionization 
cross sections. One is the uncertainty in the electron degra- 
dation spectrum. The second is the uncertainty in the rate 
for an individual process. For hydrogen these are explicitly 
coupled, but this is not necessarily so for other species since 
such species have a negligible influence on the degradation 
spectrum. In the six test cases, it turns out that the model 
spectra show no sizable difference. The Ha profile shows a 
slightly stronger dependence on the impact cross sections 
of hydrogen. However, the two approximations adopted to 
compute the excitation and ionization cross sections produce 
fairly good calculations for hydrogen, which reduces our con- 
cern about such uncertainties. Generally, a factor of five in 
the uncertainties of the cross sections will not dramatically 
influence our results. 



8.2 E max for high energy electrons 

Radioactive decay from 56 Ni and 56 Co generally produces 
fast electrons with energy ~ 1 MeV. However, in the non- 
thermal solver, we assume a source function that injects all 
electrons with an energy E max = 1 keV, which is very small 
compared to reality. While this high energy cut also intro- 
duces another source of uncertainty, the Bethe approxima- 
tion, the Arnaud & Rothenflug formula, and the electron 
thermal term, have similar asymptotic behaviors at high en- 
ergy, so the uncertainty should be small as long as E max is 
large enough. We have run two additional models with E max 
= 2 keV and E max = 0.5 keV. The difference in the fractional 
energies entering the three channels between one of the test 
models and the model with E max = 1 keV is less than 5%, 
and the difference between the two test models is slightly 
larger. Importantly, the choice of E max makes no difference 
to the spectra, de spite the effects on the deposition fractions. 
Previous work bv lXu fc McCravl (|l99ll ) showed that the de- 
position fractions only have a weak sensitivity to E max . 



8.3 The time-dependence effects on non-thermal 
processes 

Another uncertainty comes from the late inclusion of non- 
thermal excitation and ionization. At very early times the 
region where 56 Ni decays is in LTE, and the assumption 
that all the energy goes into heating is excellent. However, 



we included non-thermal excitation and ionization on day 
40.6, when the time-dependent effect is already functioning 
at large velocities, in regions of low ionization and under- 
abundant in 56 Ni/ 56 Co. We study this uncertainty in two 
ways. The first method is to compare the thermal and the 
non-thermal model at day 40.6. At this epoch, non-thermal 
excitation and ionization are still unimportant. The model 
spectrum is hardly affected by the inclusion of non-thermal 
processes. This is also the main reason why we start the 
non-thermal sequence at this epoch. The second method is 
to compare two non-thermal models at a later time - one 
is evolved from day 40.6 and the other with immediate in- 
clusion of non-thermal processes. At day 127, we found very 
subtle differences between the spectra of the two models. 
The Ha profile is only slightly weaker in the model with im- 
mediate inclusion of non-thermal processes. The above two 
tests indicate that the late inclusion of non-thermal excita- 
tion and ionization only introduces very small effects on the 
nebular spectra. Future studies will start with non-thermal 
processes, thereby removing such uncertainty. 



9 DISCUSSION 

The amount of mixing i n SN 1987A has been widely debated 
(See lArnett et ai]|l989l for a review). The two fundamental 
questions are what velocity hydrogen is mixed down to, and 
what velocity 56 Ni i s mixed outward to. T he latest 3D sim- 
ulations of CCSNe (|Hammer et al.ll2010t ) showed that sig- 
nificant amounts of hydrogen are transported into deep lay- 
ers of the ejecta 1000 kms -1 ), and that a large amount 
of heavy elements are strongly mixed above 2000 kms - . 
These results are thought to be very similar to the case of 
SN 1987A. However, it is difficult to draw a conclusion, since 
the simulations hav e a diffe rent progenitor from SN 1987A. 
iKozma fc Franssonl l|l998bl ) modeled the late time spectra 
of SN 1987A and found hydrogen was mixed down to ^ 
700 kms" 1 (they needed th e mixing to fit the observations) . 

iFassia fe Meikld (|l999h modeled the He 1 1.083 /im line 
and found some amount of 56 Ni is required to be mixed 
above 3000 kms -1 to obtain a best fit. lMitchell et al.l (|200lh 
synthesized the spectrum of SN 1987A and argued that the 
strong Balmer lines at early times is due to non-thermal ef- 
fects, which required a substantial amount of 56 Ni mixing 
out to 10 000 kms -1 in the hydrogen envelope - velocities 
much lar ger than predicted i n hyd rodyn amic simulations. 
Howev er, Tutrobin fc Chugail (|2005l 'l and iDessart fc Hillierl 
l|2008h showed time-dependent ionization is the main reason 
for the strong Balmer lines. The inward mixing of hydrogen 
and outward mixing of heavy elements in SN 1987A are still 
open questions. 

The scale of mixing also plays an important role. Mi- 
croscopic mixing results in homogeneous compositions in the 
ejecta, while macroscopic mixing can lead to large-scale com- 
positional inhomogeneities. With microscopic mixing the 
mixed species can interact directly (e.g., charge exchange 
of H with Fe) . With macroscopic mixing, this will not occur 
even if H and Fe are at the same radius (but different 8, <f>). 
Microscopic mixing can be treated in ID, while an accurate 
treatment of macroscopic mixing may require 3D models. In 
SN 1987A, the high velocity feature of ir on-group emissio n 
lines (e.g., [Fen] 17.94 ^m and 25.99 /im (|Haas et al.lll990l '); 
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[Nill]6.6/xm jColgan et al.l Il994)), at late times s t rongly 
indicate macroscopic mixing. iFransson fc Chevalier] 0989) 
showed that microscopic mixing can significantly ch ange the 
observed spectra at late times, and lLi et al.l (|l993T l empha- 
sized the importance of a clum py structure on the Fe/Co/Ni 
lines. Ijerkstrand et all (|201ll ) were able to reproduce rea- 
sonably the SN 1987A spectra 8 years after explosion with 
a purely macroscopic mixing model. Furthermore, hydro- 
dynamical mod els were unable to produce e fficient micro- 
scopic mixing (fshigevama fc Nomotol ll99(J ; iMueller et alj 
ll99ll : lHerant fc WooslevHl994l ). 

Modeling the late time spectra of SN 1987A provides 
insight into the amount of mixing. The disappearance of 
Balmer continuum photons at the nebular phase makes 
non-thermal excitation and ionization the only way to re- 
produce the Balmer lines. However, the problem is com- 
plicated by the 7-ray transport. 7-rays interact with the 
medium through three processes: photoelectric ionization, 
Compton scattering and pair production. The decay of 56 Ni 
and 56 Co mainly produces 7-ray with energy of ~ 1 MeV. 
At this energy range, Compton scattering is the dominant 
process. To give a rough estimate of the mean free path of 
a high energy electron before it is Compton scattered by 
atoms, we adopt an effecti ve, purely absorpti ve, gray opac- 
ity k 7 = 0.06Ye cm 2 g _1 jSwartz et all 1 19951 ). where Y e is 
the total number of electrons per baryon. The optical depth 
is calculated by, 



dr — —K-ypds 



(14) 



Assuming Y e — 0.5 and integrating r from velocity 
1500 kms" 1 , the place where 56 Co is still abundant, we 
obtain the mean free path of 7-ray is ~ 1.0 x 10 14 cm, 
which means the photon can travel from velocity 1500 to 
1700 kms -1 in model D127.NT. No doubt, very few 7-rays 
travel a long way before they interact with the medium and 
deposit their energies. A Monte Carlo co de has been devel- 
oped for computing the 7-ray transport (|Hillier fc Dessart] 
I2012T ), which will give a more detailed understanding of the 
mixing problem. We ran a test with allowance for 7-ray 
transport at day 127 and, as expected, there is only a slight 
enhancement in the line strength and width of Ha. 

The first discovery of ha rd X-rays in SN 1987A was 
made bv lSunvaev et all (|l987i ). approximately 5 months af- 
ter th*_e2cpl^sioii J _and the first detection of 7-ray was made 
bv lMatz et al.1 l)l988h . shortly after the discovery of the hard 
X-rays. This suggests that 7-ray transport was beginning to 
become important at day 127. However, our model at day 
154 gives an estimate that 7-rays can only travel from 1500 
to 1800 kms -1 using the effective gray opacity re 7 , which 
indicates an insufficient outward mixing of 56 Ni in the in- 
put. In fact, theoretical modeling of the bolometric, X-ray 
and 7-ray light curves require subst antial outward mixing 
of 56 Ni into the h ydrogen envelope (|Kumagai et al.l 1 19881 ; 
lArnett fc FuHl989T l. 



10 CONCLUSION 

We developed a non-thermal solver, which takes into account 
ionization and excitation due to non-thermal electrons cre- 
ated by 7-rays that arise from the decay of 56 Ni and its 



daughter products, and incorporated it into our fully non- 
LTE time-dependent code to model SNe. We use the model 
l lml8a7Ad' of Woosley on 0.27 day, enforced into homology 
as an initial input, and benchmark the non-thermal solver by 
comparing model results with the observations of SN 1987A. 

Non-thermal models have lower temperature and more 
excited/ionized material in the region where the non- 
thermal processes are crucial. Fractional energy of non- 
thermal excitation/ionization is prominent at places with 
very small ejecta ionization fraction X e . The spectral com- 
parison with the thermal models shows that Ha is strongly 
increased at nebular epochs. At late times, He I lines, which 
are totally absent in the thermal models, are present in non- 
thermal models. While most He I lines are significantly con- 
taminated by other lines, He 1 2.058 pm provides an excellent 
opportunity to infer the influence of non-thermal processes 
on helium. Hel7065A is a possible optical line that can be 
used to infer the influence of non-thermal effects in the nebu- 
lar epoch. There are many Ol lines locating near 1.129 /im, 
and they are also strengthened by non-thermal processes. 
Most of other lines are only weakly affected at the epochs 
considered here. 

Optical and IR Hel lines mainly originate below 
2000 kms" 1 . We confirm that non-thermal excitation is the 
most important process for Hel 2.058 fim. However, Hel 
1.083 /jm is due to cascades from higher levels, which indi- 
rectly relate to non-thermal ionization. Although photoion- 
ization and recombination are prominent processes for pop- 
ulating many levels, non-thermal excitation and ionization 
are the processes controlling the ionization balance. 

We also compare the non-thermal models with the ob- 
servations. The re-brightening of the bolometric light curve 
peaks about 10 days later in our models. An underestimate 
in the amount of outward mixing of 56 Ni is the mam reason 
for the delay. Although the synthetic bolometric light curve 
agrees reasonably with the observations, multi-band light 
curves show discrepancies in various degrees and reveal po- 
tential differences between the hydrodynamical model and 
SN 1987A or possible problems in our modeling. The opacity 
issue, related to limited levels in the atomic data or some 
missing species, is one possibility, considering that the V- 
band is overestimated and the I-band is underestimated. 
Ha maintains a strong profile at the nebular epochs in the 
non-thermal models, resembling observations of SN 1987A. 
However, the model Ha profile is much narrower than what 
we observed in SN 1987A, probably due to insufficient out- 
ward mixing of 56 Ni. Although hydrogen is abundant above 
2000 kms -1 in the model, the lack of 56 Ni in this region 
and the assumption of local deposition of 7-rays make the 
non-thermal effects negligible there. We also compared the 
optical spectral evolution to the observations and found a 
satisfactory agreement. 

The uncertainties introduced by various approximations 
and assumptions in the non-thermal solver have been in- 
vestigated. The primary uncertainties are due to the not- 
well-known electron impact excitation and ionization cross 
sections of metal species. However, our results for a BSC 
model are largely dependent on accuracies of the cross sec- 
tions for hydrogen and helium, since they are the most 
abundant two element s. The Bethe approximation and the 
lArnaud fc Rothenflud formula produce very good estimates 
of the excitation and ionization cross sections for these two 
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elements. Therefore, our conclusions are not affected by 
these uncertainties. Another source of uncertainty comes 
from the computation of the degradation spectrum, which 
arises from the choice of the number of energy bins and the 
energy cutoff E max at the high end of the degradation spec- 
trum. We demonstrate that our choice for the number of 
energy bins, N = 1000, is sufficient to give a reliable degra- 
dation spectrum and using E max = 1000 eV in our model 
produces almost the same spectra as using E max = 2000 eV. 
We also explore the time dependence of non-thermal pro- 
cesses and find there is only a very subtle effect on predicted 
spectra. 

With the non-thermal solver, we are able to simulate 
SNe from photospheric to nebular phases continuously. The 
influence of non-thermal ionization and excit ation on Type 
lb an d Type Ic models is being investigated jDessart et al.l 
2012). As more and more nebular spectra are available, the 
comparison of observations with models will allow us to 
place constraints on the hydrodynamic models and nucle- 
osynthesis. The non-thermal solver also provides opportuni- 
ties to constrain mixing effects in SNe. 
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